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ABSTRACT 


The development of a numerical procedure for the treacke 
ment of nonsimilar, unsteady, laminar boundary layers is 
presented. The solution is obtained from the laminar, iso- 
thermal, incompressible boundary layer equations employing 
a modification of the integral matrix procedure of Bartlett 
and Kendall. Solutions of example problems are presented 
for steady Blasius and Howarth flows, and for oscillating 
Blasius flow. Agreement with the known classical results is 
satisfactory and establishes the general feasibility of the 
method. Core storage requirements of 130,000 bytes allow 
consideration of as many as 25 nodal points across the 
boundary layer, 50 time increments per oscillation cycis 
and 50 streamwise stations. Solution of oscillating Blasius 
flow considering 8 nodal points and 16 time increments 
requires 13.49 seconds for one streamwise station utilizing 


IBM 360/67 time sharing capabilities. 
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t. SiNTRODUCTION 


A. GENERAL 

All aerodynamic flows are to some degree influenced by 
unsteadiness. Rotating stall in turbomachinery, aerodynamic 
stall flutter and hydrofoil flow fields are but a few 
practical problems in which the nonsteady fluid flow phe- 
nomena are of prime importance. Additionally, the recent 
emphasis on the development of rotary wing aircraft has 
created a new demand for knowledge of the behavior of 


boundary layers in unsteady flow. 


B. ANALYTICAL HISTORY 

Early analytical treatments by Stokes, Rayleigh and 
Schlichting (Ref. 1) considered only a special case of the 
general problem: that of unsteady viscous fluid flow with 
no mean flow or pressure gradient. Lighthill (Ref. 2) made 
the first significant investigation of small sinusoidal 
oscillations superimposed on a mean flow. Lighthill's 
analysis was based on a small perturbation treatment of the 
boundary layer equations for sinusoidal flow, with only 
the first order terms retained. The determination of the 
existance of the "quasi-steady" regime at low frequencies, 
and the "shear wave" regime at higher frequencies was the 
principal result of this analysis. Glauert (Ref. 3) extended 
Lighthill's work by considering the boundary layer in the 


vicinity of a stagnation point. Carrier and DiPrima (Ref. 4), 


dpe: 


using a system of equations derived using a modified Oseen 
linearization, confirmed the existance of the "Shear wave" 
solution and the phase advance of the wall shear determined 
by Lighthill. Nickerson (Ref. 5) initially retained per- 
turbation terms beyond the first order, thus widening the 
scope “Of Lightehi Mesework. 

In addition to his analytical development, Nickerson 
made hot-wire measurements of the laminar boundary layer on 
a flat plate oscillating in a Blasius mean flow. His work 
partially confirmed both his own and Lighthill's analyses, 
but mechanical difficulties precluded the completion of the 
experimental program. Hill and Stenning (Ref. 6) imposed 
Sinusoidal oscillations directly on mean Blasius and Howarth 
flows and were able to confirm the analytical results of 
both Lighthill and Nickerson. Recently an integral analysis 
of the unsteady boundary layer was introduced by Koob and 
Abbott (Ref. 7), who succeeded in demonstrating the method 


for the case of an accelerating flat plate. 


C. APPLICATION OF NUMERICAL METHODS 

With the advent of the high-speed digital computer, it 
has become possible to apply numerical methods to the 
analytical studies of viscous fluid flows. Fromm (Ref. 8) 
and Chorin (Ref. 9) introduced modified finite difference 
methods for the treatment of viscous flows which could 
feasibly be applied to the boundary layer problem. Smith 


and Clutter (Ref. 10) and Kleinstein and Nabi (Refs. 11 and 
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12) applied finite difference methods to the steady boundary 
layer equations. Successive approximation, linearization 
and integral methods have also been applied to the boundary 
layer equations (Refs. 13, 14 and 15). A major drawback of 
most of the methods mentioned is that in order to reach a 
solution to a given problem large computing times are nec- 
essary. 

The introduction of the Boundary Layer Integral Matrix 
Procedure by Bartlett and Kendall (Ref. 16) provided a major 
breakthrough in the field of fast computer solution of the 
steady flow boundary layer problem. The method was developed 
as a means of obtaining rapid computations in the laminar 
regime for steady, multicomponent, chemically reacting 
boundary layers; and the logical extension of the technique 
was to encompass the unsteady flow condition. 

The purpose of this investigation is to develop a com- 
puter program capable of solving the unsteady, laminar 
boundary layer problem by adapting the rapid solution tech- 


niques of the Boundary Layer Integral Matrix Procedure. 
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If. ANALYSIS 


GOVERNING EQUATIONS 


In this development, utilizing the integral matrix 


method, it is convenient to make the assumption of incom- 


pressible isothermal flow. The resulting laminar boundary 


layer equations for two-dimensional (n=0) or axisymmetric 


(n=1) flow are as follows: 


Continuity: 
n 
one) 

Xr = y 

oO 
Momentum: 

gu, du, pul _L oP yu du Ss 

ot os oy 0 os 0 a2 

aU CUT. ee Oe 

Saeed fh) 


where s and y are respectively the streamwise and normal 


components, u and v are the corresponding velocity compo- 


nents, U is the streamwise velocity component outside the 


boundary layer, ry is the radius of the body perpendicular 


to the axis of revolution for an axisymmetric body, n is 


zero for a two-dimensional body and unity for an axisymmetric 


body, t is time, p is density, uw is viscosity and P is 


pressure. 
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B. BOUNDARY CONDITIONS 
The boundary conditions considered for the governing 


equations as stated are: 


Mpen te 0 (4) 
Mies 0 (5) 
Un U's ; &) (6) 
dU, 


where the subscripts w and e refer respectively to wall and 
boundary layer edge conditions. Equation (4) is the require 
ment of zero slip at the wall, equation (5) implies no mass 
transfer at the wall, and equations (6) and (7) are repre- 
sentative of the assumption that the boundary layer edge 
velocity is equal to the free stream velocity. 


The free stream velocity may be written: 


U = (U, = bys) (1 + € cos wt) (8) 


where Un is the mean free stream velocity at the leading 
edge, by is the rate of change of the mean free stream 
velocity with distance from the leading edge, € is the ratio 
of the amplitude of the free stream oscillation to the local 
mean free stream velocity and w is the frequency of the 


oscillation. 


ie 


C. TRANSFORMATION OF EQUATIONS 

Paralleling the development of Bartlett and Kendall 
(Ref. 16) a transformation of variables is performed util- 
izing a combination of the Levy and Mangler and the Howarth- 
Dorodnitsyn transformations which is known by a variety of 
names including Levy-Lees, Lees-Dorodnitsyn and Mangler- 


Dorodnitsyn. This transformation is as follows: 





: 2n 
—E=f UsPauat es ds (9) 
O 
Ce y 
n = i peasy (10) 
V 2° GeO 


Simplifying equations (9) and (10) to encompass the 
assumptions of incompressibility and constant temperature 


the transformations become: 


= 2n 
Exv= Jour Wee ds (15) 
oO 
n 
Foe 
¥Ze 


The transformation to Levy-Lees variables is greatly 


simplified by use of the operator: 


of 


oar on 








Q 
— 
alll 
ww 
~_- 
Fe 
WwW 
all 


oS vet} = atu_r” /2E aS 
S eo 


( (ee 
oS 0& Ze 


where £ is the stream function defined as: 
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f-f =f —dn (14) 
O 


the prime represents partial differentiation with respect to 


n so that: 


 -» 7 
gh ae (15) 
Se 


and a* is a normalizing parameter: 


n 
pyu_r 

os re (16) 
V2E 


In addition, the partial derivatives are: 











GM, eat a(_) , an ac) 
ae = a*r Y2E { OE os Oe on } (17) 
je. aes | 
we SO) (18) 
2 2 ww 

2) = (SE y" (19) 
oy 


Applying equation (15) to the velocity in the partial 


derivative with respect to time yields: 


§ 
aie 4 2 (20) 


~ Ue. Ze at 


cle 


Using equations (13) through (20) and applying equation 
(6) to equation (3) results in the following transformed 


momentum equation: 


a9 


tet i ee 
ft V£f- tte =o 








' : d2nu | 
a ot ee BOE) Sere iets oy ey,_ | 
21f sme 7! ine’) =o (oe | Se eee 
' 
where | 
ginu, 


Transforming the boundary conditions, equations (4) 


through (7) become: 


fa (23) 
ae (24) 
te (25) 
fa (26) 


D. TAYLOR SERIES EXPANSIONS 

At any streamwise station s, the boundary layer may be 
divided into N-l strips which are joined by N nodal points 
Nye where iis 1 at the wall and N at the boundary layer 
edge. The primary dependent variables f; and their normal 
derivatives are related by Taylor series expansions such 
that the eh are represented by connected cubics with con- 
tinuous first and second derivatives at the nodal points. 
This is commonly called a spline fit. 


Expanding the primary variables in Taylor series and 


tre 
truncating rat f. the following linear equations result: 
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oO pee Se GT 


-f£,,,+£,+#£, n+ f£, Dee, sit = 0 (28) 
SEA ee tee e = yey ~ = 0 (29) 
oes Ve one (30) 


Equations (27) through (29) represent a set of 3(N-1) 
equations where f is represented as a quartic, f as a cubic, 
f' as a quadratic and f£'' is considered to vary linearly 


between ns and Nesd° 


E. AXIAL AND TIME DERIVATIVES 


The axial derivatives of equation (21) are treated as 


logarithmic two or three-point differences such that: 








fac) _ 
eerie = oe 1 et y= heey 
is a 
where for two-point differences: 
i= d, = -d d, = 0 (ep 
O may i O 2 
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and for three-point differences: 


gag-1te4e-2 grg-2 gog-1 


d= 2 OE" 22> 


Ophea 222-2 + “ghoey gaits  § 2 yee 
(34) 
in equations (32) through (34), ( ) gjimplies the previous 


streamwise station. 


A 


Differencing the time derivatives in the same way: 


Penh = 
eopegl eee De A ey oo 


where for two-point differences: 


T = f= \-= i T= 0 (37) 





and for three-point differences: 


_ kona tk k-2 ____k4x-2 ___k’k-1 


Lu T. = 
O A A 7 A 2 caK-2 k-1 K-2 


Z q 
Ko x-1 Ko K-2 eo k=l el poe 


(arch 
In equations: (26) through (38), ( )-jimplies the previous 


time. 


A =e (39) 


Ze 


Equations (37) through (39) may be greatly simplified 
by considering the utilization of equal time spacings. The 
resulting equations for the treatment of time derivatives 
over equal time spacings are as follows: 

Two point differences: 


- a2 = 
T= ne T) = To T. 0 (40) 


Three-point differences: 


2 : 


gs i ee an 
aes oe i aoe DS ae (41) 
where: 
eee ee eee (42) 


F. INTEGRATION OF THE BOUNDARY LAYER EQUATION 
The transformed streamwise momentum equation, equation 


(21), is integrated at constant € from eS £o Nas yielding: 


ul 


pap i ‘ a 2 
| =] + f ff dn + Bf (1 -="£ ) dn 


i-l i-l i-l 











1 Of! “OF ou 3 «ODE 
-—2 f (f ~=-f ~dn - f w= dn 
aT dene dene ee oe 
at i d£nu, 
= oD oe (f - 1) xe dne=-0 (43) 
a 1=1 


Expanding equation (43) in Taylor series and applying 


equation (36) to the time derivative terms: 


a3 








r of! dn a / a | 
i-l ot _ét aa 
(Tf, + Tif ,-1 + Tof, 2] (44) 
i d2nu dknu. r i 
Paige =a, dn = (ie] - én) 
eal ot ot LL are 


i. 
= [t | -6ni[T 2nu_ + T.nu +T. &nu 1(45) 
1 oe Ea ee 


the integrated boundary layer equation becomes: 


aat 


Pi 
g-1 ts eee rite +7 


+ +f (l+d jf + (d,f 0 tT} k-1 tT oF, 2! 


i 
-f£[T &nu_ + T. Snu + T. 2nu | 
Oo e e C,_9 


er t Sy) 2 =a 


-(1+ B+ 2d.) [f; XM, + f£, XM, + £; XM, + £;_)XMj,] 


-2[f. 2M. + £° 2M. + £:°2M. + £:'\ 2M 
Bf 2 a i 





all 2 3 1- 
OU 
+ + = 
6n{B + ee [T &nu, T)&nu, + T,%nu, ae 0 (46) 


The Taylor series expansions, XM and 2M, of the integrals 


in equation (43) are listed in Appendix A. 


G. RECURRENCE FORMULAS 
The boundary layer equations, equation (46), and the 
Taylor series expansions of the primary variables, equations 


(27) through (29), together with the boundary conditions, 
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equations (23) through (26), are solved by Newton-Raphson 


iteration of 4N simultaneous equations in 4N primary 





variable unknowns. These are summarized as follows: 
Primary Variables No. of Equations 
is N 
a N 
if N 
vas 
£ N 
No. of Eqns. Req: 4N 
Equation Nos. No. of Equations 
Taylor Series Expansions (27029) 2 (Net) 
Boundary Layer Equations (46) N-l 
Boundary Conditions (34)- (26) 4 
No. of Eqns. Available: 4N 


The Newton-Raphson solution technigue is illustrated by 


considering two simultaneous nonlinear equations: 
F(x,y) = 0 G(x,y) = 0 (47) 


which have the solution x=X and y=Y. Defining Xn and Yn 2s 


th 


the values of x and y for the m iteration, the desired 


solution f£(X,Y) is expressed in the Taylor series expansions: 


OPtx  ,yi) 
CRC ey eae ae CN ce) ——s5 m 
tay ange oe : 
0 = G(X,Y) = G(x, 1Y,) + (X - x1) een’ Yn’ 
+ yee) nen, nee (48) 
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In the Newton-Raphson method (X,Y) is replaced by 
(xe Yme) and the nonlinear terms in x and y appearing in 


equations (48) are ignored yielding the set of recurrence 


formulas: 

Ax, naa + AY, aoa = ~G(x_+Y,,) (49) 
where: 

AX = Xmt1 7 Xm OYia = Smee ela (50) 


The Ax. and Ay. are corrections to be added to x. andy 

™m m m m 
respectively yielding Xt] and Vint’ F(x +Y,) and G(x, +Y,) 
are the values of the original functions F(x,y) and G(x,y) 


evaluated for the value of the variables during the meh 
iteration. Noting that Po and G(x, +Y,) approach zero 
as the corrections approach zero, it is appropriate to con- 
Sider them as errors associated with the original equations 
(47). 

Differentiating the Taylor series expansions, equations 
(27) through (29), with respect to the primary variables in 
accordance with the derivation of equation (49) yields the 
following recurrence formulas for eens iteration: 

Z 5 5 
(-1) AE. .,+(1) AE. +(6n) AE) + (Sy ne + (Sy ae"+(22_) ne!" =- ERROR 
itl i at 2 al 8 a 24 a 


(oi) 
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2 2 
- ' ' " én wy 6n un io 
( L) AE; 4 + (1) AE; + (6n) AES +(——-) AE; tee) oe ERROR (52) 


foal) Acameee (1) Ae, + (Dy ae,” + (By ats = - ERROR (53) 
where Af. represents the correction for fis the coefficients 
of the corrections represent the partial derivatives of the 
Taylor series expansions with respect to the primary var- 
lables, and the ERRORs are equations (29) through (31) 
respectively, evaluated for the values of the variables 
during the iteration. 

The recurrence formulas for the boundary layer equations 
are: 


at “te /Nels (ltd )£ = f(l+d )At SAVE (djf,_4+ dof) _5) 





, aE 
__ PU = i 
ae A£(T gnu, + T)<nu, + Tj knu, ) } 


oO ee en 
' W vee uns 
-2(1+ 8 +2d.){Af,xM,+Af,xM, + Af, XM, + Af. _,XM,} 
' " we 
-2{Af, 2M, + Af, 2M, + Af, 2M, + Af,_, 2M,}= —- ERROR (54) 


where ERROR is the left-hand-side of equation (46) evaluated 
for the values of the variables during the iteration. 
Finally, the recurrence formulas for the boundary con- 


ditions, equations (23) through (26), are: 


Af = - ERROR 
W 


- (£) (55) 


Af = - ERROR = ie) (56) 


27) 


Af. = - ERROR = l - a (7) 


Af = - ERROR 


eee (58) 


Equation (54) may be further refined by collecting terms 
and writing the equation in terms of correction coefficients 


as follows: 


CLA£. + C2Af. . + C3AF: + C4AF: 
a a= 1 L= 


il a 


Ty " vit wet 
+ CSA£; + C6Af,_, + C7Af, + C8A£,_,=- ERROR (59) 


The. correction coefficients €l throvnghece are g2veneam ' 


Appendix B. 


H. MATRIX FORMULATION 
The coefficients for the recurrence formulas, equations 
(51) through (53) and (55) through (59), evaluated for the 
hi 


m iteration form a square matrix [a with 4N rows and 


columns. The matrix equation is: 
a Wea le |= | (60) 
By oll aaa 
where the column matrix [ av | is the matrix of corrections to 
the primary variables, and the column matrix [=] is the 
matrix of errors. The matrix [a] contains many zeros, and 
great savings in computer storage space and computation time 


may be realized through proper ordering of the equations and 


variables. 
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The equations are first divided into linear and nonlinear 
sets, and at the same time, the variables are classified as 
linear or nonlinear, resulting in the partitioned matrix 


equation: 


soem (61) 


werent ee 
! 
1 
i 
1 
' 
( 
‘ 
' 
( 
‘ 
‘ 
‘ 
‘ 
| 
t 


The division of variables into linear and nonlinear sets is 
generally arbitrary except that the division must be made to 
insure that the square matrix AL is not singular. Adopt- 
ing the variable groupings of Bartlett and Kendall (Ref. 16) 


the linear and nonlinear corrections are arranged as follows: 


' ' ' " tt ve 
AVL (Af, ,Af3,...,Af_,Af At poee AE, AL ,AfS,.--,Af,_)) 


eee, ) 


AVNL (Af ,Af.,Af , ar if 


The linear and nonlinear equations are sequenced within the 
partitioned square matrix EN as shown in Figures 1 and 2, 
which demonstrate the case of 4 nodes across the boundary 
layer. 

Expanding equation (61), the resulting linear and non- 
linear equations are, respectively: 


f 7 iy; We 1 
fat] vt lea ae] avn |= [en | (62) 
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as ne : 
jane | avn | | put | avwt, | = “ent | (63) 
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Solving equation (62) for the linear corrections: 
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FIGURE |.LINEAR PORTION OF [A] MATRIX 
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FIGURE 2. NONLINEAR PORTION OF [A] MATRIX 
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AVL = - AL! ~BL:; AvVNL ;}+! ani | -EL | (64) 
. » Po ee eee spo ot | eae zs 
Examination = the AL and "BL | matrices reveals that 


io te 


the product Tan? Tor] may be i ae given only the nodal 


- 


Spacing and must be determined only once for a given problen, 








thus: 
eae 
| ov | - -| Bal | over | + | axa | (65) 
L i LC 
where: 
a = aes yr 
BAL ! = j AL] La (66) 
i aul —_ a 
mS : yr 
ELA | = AL! EL, | (67) 
a se el 


Applying equation (65) to equation (63), and solving for 


the nonlinear corrections: 


/ aie i Le 
AVNL | = | Br | ENE (68) 
where: 
= a i = 
BNL «= oe ant |) aay | (69) 
= feo ae 
ENL =. Seer a NE | Re (70) 


LT. SOLUTION PROCEDURE 


The matrix solution proceeds according to the following 


algorithm: 
=1 
1. Given the nodal spacing determine [ AL | and form 
the product | BAL]. 
2. Given_initial values of the primary variables deter- 


mine "EL | and [- -ENL |, and form the product [ ELA]. 


ee a = 
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3. Form_! BNL, and _ ENT | from the coefficients of 
om] f L 

and | BNL.. 


| ANL | 
ira wl r 7 sae 2 

4. Invert | BNL, and form the product “BRE | ' ENie!, 

which is thé solution for BAe a r 


5s Selve for [ave] from equation (64) using | ave. | 
as determined in step 4. 


6. Add the linear [av] to the corresponding primary 
variables to complete the iteration. 


7. Compute the new errors '-EL | and | -enz, |, and 

test for a maximum allowable error. if this cretberien 
has not been reached, return to step 2 using these 
errors to commence the new iteration. 
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III. RESULTS AND DISCUSSION 


The computer program has been used to investigate three 
basic categories of boundary layer problems using a flat 
plate model; steady Blasius flow (zero pressure gradient), 
steady Howarth flow (uniform adverse pressure gradient), 
and oscillating Blasius flow. The velocity profiles which 
resulted from these solutions are compared to previously 


reported results and presented in Figures 3, 4, and 5. 


A. STEADY BLASIUS FLOW 

In Figure 3, computed values corresponding to three 
nodal spacing choices are superimposed on the Blasius pro- 
file. Agreement between the numerical solution and the 
classical solution is better than 0.08 percent for as few 


as 10 nodal points. 


B. STEADY HOWARTH FLOW 

Results for a steady flow in a uniform adverse pressure 
gradient are compared with the solution of Howarth (Ref ig) 
in Figure 4. The comparison of the computer solution with 
that of Howarth is accomplished by correlating the ainetem 


Sionless parameter x* (=b,s/U,). which characterizes 


Howarth's velocity profiles. Agreement between the two 
solutions is generally good with differences no greater than 
4.35 percent for 22 nodal points. The differences encoun- 


tered may be reduced somewhat by increasing the free stream 


velocity of the test problem. This results in larger 
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| _ BLASIUS SOLUTION (REF.1) 

ZL. _UBLIMP~1IO PT. SOLUTION 
© UBLIMP- I2 PT SOLUTION 
O 


UBLIMP -14 PT. SOLUTION 
50 = 


40 
3.0 


2.0 





VELOCITY RATIO —(f'= U/ug) 


FIGURE 3. STEADY BLASIUS FLOW DATA COMPARISON 
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stream-wise pressure differences in the equations of motion 


and consequently lesser roundoff errors. 


C. OSCILLATING BLASIUS FLOW 

The unsteady capabilities of the numerical procedure 
were investigated in the solution of Blasius flow with 
Sinusoidal oscillations superimposed on the mean flow. 
UEttazing the critical oscillation frequency We (=.6U/s, 
for a flat plate) of Lighthill's analysis (Ref. 2) to 
determine the solution approximation with which to compare 
the numerical results, the profile of velocity fluctuations 
which resulted was compared with the "quasi-steady" 
analysis as shown in Figure 5. Agreement between the calcu- 
lated velocities and the Lighthill profile is within 15 
percent deep in the boundary layer and within 5.5 percent 
in the outer regions. Noting that Lighthill's analysis is 
based on linearized small perturbation theory and that the 
amplitude of oscillation in the test problem of 0.1 Ne can- 
not be considered a small perturbation, better agreement 


could scarcely be expected. 


D. NUMERICAL PROCEDURE 

The computer program has been operated on both IBM 360/67 
and CDC 6600 computer installations. Core storage require- 
ments are 130,000 bytes for a maximum array size encompassing 
25 nodal points across the boundary layer, 50 streamwise 


stations and 50 time increments per oscillation cycle. A 


Sy 


——— LIGHTHILL (REF. 2) 


@= 
va. WO HZ. 
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RATIO OF VELOCITY FLUCTUATIONS (Au, /AU,) 


FIGURE 5. OSCILLATING BLASIUS FLOW 
DATA COMPARISON 
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typical 25 nodal point Blasius solution was performed in 
7 seconds on the CDC 6600. 

Limitations in the present numerical capabilities of 
the program precluded a more extensive investigation of 
parametric effects, however the results qualitatively 
demonstrated the feasibility of the computational procedure. 
Further refinements to the computer program should enable 
a more extensive investigation to include quantitative com- 


parisons with experimental results. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


The feasibility of the Unsteady Boundary Layer Integral 


Matrix Procedure in the treatment of simple steady and 


unsteady flow problems has been demonstrated. Favorable 


comparison with classical results provides a firm base for 


continued investigation in the unsteady regime. 


the 


Continued research utilizing the program should include 
following revisions: 


1. Alteration of the present time derivative component 
arrays (FTM1) and FTM2) to insure that all values of the 
components are stored for each nodal point, instead of 
retaining only the values at the wall. 


2. Incorporation of Blasius profiles as first estimates 
for the primary variables in order to decrease the 

time required to complete the first nodal iteration 
procedure. 


3. A study of array usage aimed at reducing the re- 
quired storage and subsequent conversion to double 
precision computation in order to minimize round-off 
interference. 


4, Addition of the capability of handling compressible 
flows in order that higher velocities may be introduced. 


5. Inclusion of turbulent boundary layer capabilities 
perhaps through the use of eddy transport properties. 
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TAYLOR SERIES EXPANSIONS OF INTEGRALS (Ref. 16) 
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APPENDIX B 


BOUNDARY LAYER EQUATION CORRECTION COEFFICIENTS 
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43 


APPENDIX C 


PROGRAM DESCRIPTION 


A. GENERAL 

The computer program referred to as UBLIMP, is written 
in FORTRAN IV source language and has been operated on both 
IBM 360/67 and CDC 6600 computer installations. Core 
storage requirements are 130,000 bytes for a maximum array 
size encompassing 25 nodal points across the boundary layer, 
50 streamwise stations and 50 time increments per oscillation 
cycle. 

The program consists of a main program (UBLIMP) and 21 
subroutines which are divided into two sections as shown 
in Figure 6. LINK 0, common to both sections, contains the 
main program and certain service routines which are common 
to both sections in the fully expanded program. LINK 1 
sets up the boundary conditions and LINK 2 iterates to solve 
the boundary layer equations. 

The program may be applied to steady or unsteady flow, 
two-dimensional cartesian or axisymmetric bodies with similar 
or nonsimilar profiles. It allows quadratic or cubic curve 
fits of the primary variables. Program options are con- 
trolled by the control variable KR and are detailed in the 
input instructions for the UBLIMP program in Appendix D. 

A program linkage schematic is shown in Figure 7 and a 


functional flow chart is shown in Figure 8. Correlation 
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CASE 


A READ INPUT DATA FOR NEW CASE WHICH DIFFERS FROM THAT 
OF PREVIOUS CASE. 


B IF NODAL SPACING DIFFERENT FOR THIS CASE SET UP MATRICES 


FOR TAYLOR SERIES EXPANSIONS AND LINEAR BOUNDARY CON- 
DITIONS FROM NODAL POINT DISTRIBUTION. INVERT MATRIX 
OF COEFFICIENTS OF LINEAR CORRECTIONS (THIS, IN EFFECT, 
EXPRESSES LINEAR CORRECTIONS IN TERMS OF NONLINEAR 
CORRECTIONS ). 





C COMPUTE OR READ FIRST GUESSES FOR PRIMARY VARIABLES 
OR USE PREVIOUS VALUES. 





TIME 


E COMPUTE NONSIMILAR TERMS (THOSE INVOLVING STREAM- 
WISE DERIVATIVES). 


ITERATION 














F EVALUATE LINEAR ERRORS. 


G EVALUATE ERRORS AND COEFFICIENTS FOR NONLINEAR 
EQUATIONS AND CORRECT NONLINEAR ERRORS FOR 
LINEAR ERRORS. 


H INVERT MATRIX OF NONLINEAR CORRECTION COEFFICIENTS. 


| EVALUATE NONLINEAR AND LINEAR CORRECTIONS. 
CORRECT PRIMARY VARIABLES AND TEST ERRORS FOR 
CONVERGENCE. 


J PRINT OUTPUT DATA FOR CURRENT STATION AND TIME. | 
K INITIALIZE TIME. | 


STOP 


FIGURE 8. FUNCTIONAL FLOW CHART OF UBLIMP PROGRAM 
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between the two is indicated by letters running from A 


EAroughnis « 


B. PROGRAM FUNCTIONS 
1. Main Program (UBLIMP) 

The main program provides linkage between the 
boundary conditions section and iterative solution section 
of the program. Initial core zeroing and the call for 
solution output are also initiated by the routine. 

2. Subroutine ZEROIT 

ZEROIT performs the initial zeroing of all named 
common blocks to eliminate the possibility of introducing 
random numbers in the solution. 

3. Subroutine SETUP 

SETUP is the control routine for setting up the 
boundary layer edge conditions and the streamwise derivatives 
for a new time, station or case. 

4, Subroutine RECASE 

RECASE is the control routine for the input of 
boundary layer data and initializes control variables for 
the treatment of surface discontinuities. 

5. Subroutine INPUT 

INPUT reads all the data for a problem into the 
computer. 

6. Subroutine STDETA 

STDETA permits simplification of the input data deck 
by providing standard nodal point distributions should the 


user not care to provide one. 
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7. Subroutine LINMAT 
LINMAT sets up the matrices for the Taylor series 
expansions and linear boundary conditions from the eta 
spacing and solves to express the linear corrections in 
terms of the nonlinear corrections. 
8. Subroutine MATONE 
MATONE performs operations on a column of a matrix 
of coefficients [ Bx | or on a column ©f a Matrix of errors 


| -EL | (designated E in the subroutine call list) such as to 
-1 


form [ aL | | E |. 
9. Subroutine FIRSTG 
FIRSTG computes first estimates of the primary 
variables based on the nodal point distribution should 
the user not care to provide them. 
10. Subroutine REFCON 
REFCON calculates the boundary layer edge conditions 
and sets up the wall boundary conditions. 
11.) “Subroutine SLOPO 
SLOPQ defines a cubic equation for a set of points, 
calculates the average slope at each point and integrates 
the cubic equation between each pair of points. 
22) subroutine HISTxs 
HISTXI computes terms involving the axial deriva- 
tives and time derivatives and stores those upstream quan- 


tities needed for these difference relations. 
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13. Subroutine TAYLOR 
TAYLOR calculates the coefficients in the Taylor 
series expansions of the integrals which appear in the 
integral form of the boundary layer equations. 
14. Subroutine ITERAT 
ITERAT is the control program for performing the 
boundary layer iteration. 
15. Subroutine LINCER 
LINCER evaluates errors for the linear equations 
and determines the maximum linear error. 
16. Subroutine ABMAX 
ABMAX searches a given array for the entry with 
the maximum absolute value. 
17. Subroutine NONCER 
NONCER performs that portion of a boundary layer 
iteration having to do with the solution of the nonlinear 
corrections of the nonlinear equations, computes the damp- 
ing factor and applies it to the corrections and then 
applies the damped corrections to the primary variables. 
18. Subroutine IMONE 
IMONE evaluates the coefficients of (i-1l) cor- 


rections for the jth 


nonlinear equations, where i is the 
ee nodal point in the boundary layer. 
19. Subroutine AMSET 
AMSET calculates the contributions of the [ ant. | and 


applies them to | BNE | and | ENE |. 
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20. Subroutine IONLY 


IONLY evaluates the coefficients of the pth cor- 


rections for the eo nonlinear equations, where i is the 
ee nodal point in the boundary layer. 


21. Subroutine RERAY 
= 
RERAY replaces a matrix Ic | with its inverse and 
Due 


4 
to ot 


forms the product of the inverse with another matrix 
22. SUbroutine OUTPUT 
OUTPUT prints a standard boundary layer output block 
temee: Conyerced Solution, or at the end of each iteration 


if desired. 
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APPENDIX D 


INPUT INSTRUCTIONS FOR THE UBLIMP PROGRAM 


GROUP 1 ‘CONTROL: CARD, “TITLE 


CARD 1 FORMAT (2011,15A4) 


FIELD 1 (Columns 1-20) This is the variable KR which is 
used to control the various program options. 


COLUMN 1 Determines whether a new set of ETA values 
is to be input for the present case. 


0 Uses resident values from previous case. 
1 Values input by user (mandatory for first case). 


COLUMN 2 Designates the type of first guesses to be 
utilized for primary variables. 


0 Uses built-in relations to calculate first 
guesses. 


1 First guesses input by user. 


2 Uses resident values from previous case (cannot 
be used for first case). 


COLUMN 3 Determines the treatment of streamwise 
derivatives. 


0 Performs similar solution at each streamwise 
station. 


1 Considers two-point difference relations at all 
stations except that a similar solution is per- 
formed at the first station for non-blunt bodies 
and at the first two stations for blunt bodies. 


2 Considers three-point difference relations at all 
stations except that a similar solution is per- 
formed at the first station and a two-point solu- 
tion is performed at the second station for 
non-blunt bodies; similar solutions are performed 
at the first and second stations and a two-point 
solution is performed at the third station for 
blunt bodies; and a two-point solution is per- 
formed for the first station after a discontinuity. 
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COLUMN 4 Determines when output block is to be 
printed. 


O Output block printed for converged solution or 
for nonconverged solution after 50 iterations. 


1 Output block printed after each iteration. 
COLUMN 5 Determines treatment of Entropy Layer. 
(Not used in present program. KR(5)=4 
is input to maintain input data sequencing.) 
COLUMN 6 Designates body shape. 
O Axisymmetric blunt body. 
ine P Fanar blunt body 
2 Axisymmetric sharp body. 


3 Planar sharp body. 


4 Axisymmetric or planar shape which has no sharp 
tip or blunt stagnation point, such as a nozzle. 


COLUMN 7 Steady, zero pressure-gradient noise filter. 
O Normal computations of DUDS, 
LeeDUDS is forced to zero. 
COLUMN 8 Designates form in which wall mass fluxes 
are input. KR(8) is not utilized if wall 
mass fluxes are not input. 


0 Wall fluxes input in LBS/SEC FT**2 


l Wall fluxes input in normalized form (divided by 
-ALPHASTAR) . 


COLUMN 9 Together with KR(1l1) this designates the 
type of wall boundary conditions. 


O Assigned stream function at the wall. 


1 Assigned mass flux at the wall. KR(9)=1 is used 
if zero mass flux is used. 


COLUMN 10 Determines the type of curve fits employed 
to represent the primary variables. 
(KR(10)=0 is recommended for accuracy for 
most problems, however for severe problems 
KR(10)=1 is better since the cubics can 
become poorly behaved.) 
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QO Utilizes connected cubics. 


1 Utilizes connected quadratics except for the 
outermost segment where connected cubics are 


used. 


COLUMN 11 


COLUMN 12 


Together with KR(9) this designates the 
type of wall boundary conditions. (Presently 
KR(11)=0 is the only active value and 
requires an assignment of wall temperature.) 


Provides access to standard nodal-point 
distributions built into the program. 


O User inputs data for nodal-points. 


1 8 points with more concentration at the wall. 


2 8 points equally spaced. 


3 10 points equally spaced. 


4 12 points equally spaced. 


5 14 points equally spaced. 


6 15 points with more concentration at the wall. 


7 18 points equally spaced. 


8 22 points equally spaced. 


9 25 points equally spaced. 


COLUMN 13 


Permits the assignment of a convergence 
damping factor. (This factor is over- 
ridden if a smaller damping factor is 
computed internally by some constraint.) 


0 No damping factor is assigned. 


J If J is greater than zero, corrections are 
damped uniformly by J/10. 


COLUMN 14 


COLUMN 15 


Non-zero entry causes a complete set of 
primary variables to be output for future 
use as first guess inputs. 


Non-zero entry provides debug output for 
first guesses and linear matrices. 
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0 No debug. 
1 First guesses are dumped. 


2 Linear matrices before and after inversion are 
also dumped. 


COLUMN 16 Non-zero input indicates the pressure pro- 
file has been input in the form of pressure 
coefficients. 


COLUMN 17 Non-zero entry provides debug output for 
coefficients. 


0 No debug. 


Y For (Y+1-ITS) greater than zero, where ITS is the 
number of the current boundary layer iteration, 
the coefficients which combine to make up the 
nonlinear equations are dumped and the derivatives 
of the nonlinear equations with respect to the 
nonlinear variables are dumped before and after 
inversion. 


COLUMN 18 Not used. 


COLUMN 19 Non-zero entry provides debug output for 
linear and nonlinear errors. 


COLUMN 20 Determines which of the two internal sets 
of first guesses are to be used if the user 
does not provide them. 


0 Straight line Blasius inputs are calculated for 
the stream function and velocity ratio, and zero 
is input for the first and second derivatives of 
the velocity ratio. 


1 Straight line Blasius inputs are calculated for 
all the primary variables. 


FIELD 2 (Columns 21-80) CASE. Title of the case 
(ALPHA-NUMERIC). Used for identification of 
PEinted Gulcpuc. 
GROUP 2 TIMES AND STATIONS 
CARD 1 FORMAT (I3) 
FIELD 1 (Columns 1-3, right-justified) NITEM Number of 


time grid points per cycle when considering an 
unsteady solution, otherwise input 1. (max=50) 


SS 


CARD 2 FORMAT (E10. 4) 


FIELD 1 (Columns 1-10) TIME(1) Value to be used in 
thesoutput bleck to identify a time of solutiom 


CARD 3 FORMAT (I3) 


FIELD 1 (Columns 1-3, right-justified) NS. Number of 
streamwise stations. (max=50) 


CARD 4 FORMAT (8E10.4) 


FIELD 1,2,...8 per card(10 columns per field) S(IS) 
Streamwise distance upon which the boundary 
layer solution is based in feet. Blunt body 
problems should start with S=0.0. Sharp body 
or nozzle problems must start with some finite 
value of S. The boundary layer is assumed 
similar up to and including this first station. 
A negative entry for S indicates a discon- 
tinuity at that station, and alters the differ- 
encing scheme for axial derivatives. 


GROUP 3 NODAL DATA (Skip this GROUP for KR(12)# 0) 


CARD 1 FORMAT (I3) 


FIELD 1 (Columns 1-3, right-justified) NETA Number of 
nodal points across the boundary layer. 
(maximum of 25) 


CARD 2 FORMAT (8E10.4) 


FIELD 1,2,...8 per card(10\ columns per field) ETACGE 
ETA stations across the boundary layer starting 
at the wall (ETA=0.0). It is recommended that 
a value between 5.0 and 6.0 be assigned at the 
boundary layer edge. UBLIMP values of ETA are 
equivalent to Blasius values divided by the 
square root of two. 


GROUP 4 BODY SHAPE DATA 


CARD 1 FORMAT(2E10.4) USED ONLY IF KR(6)=0 or l. 


FIELD 1 (Columns 1-10) CONE. Cone half angle in sphere- 
cone shape bodies. Leave blank for other body 
shapes. 


FIELD: 2 (Columns 11-20) RNOSE. Effective nose radius 


in feet. Used to calculate stagnation point 
velocity gradient from Newtonian relations. 
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PO a Pal 


i as 


Canwe2 FOr emil0.4) USED ONEY IF KRiG)=0, 2 or 4 


BGERLD ),2,...6 per card({(10 columns per field) ROKAP (IS) 
This is the local body radius in feet normal 
to the body centerline raised to the n power, 
where n is unity for axisymmetric bodies and 
zero for planar bodies. Therefore, ROKAP is 
unity for planar bodies and the local body 
radius for axisymmetric bodies. For planar 
bodies, this card set is used only if KR(6)=4. 
(A special input format can be used for spher- 
ical-nosed bodies as follows. Set ROKAP (1) 
equal to minus the nose radius. The nose radius 
is then set to -ROKAP(1) and ROKAP(1) is set to 
zero. If subsequent ROKAP are input as zeros, 
the program computes ROKAP from S for a spher- 
ical nose. The first non-zero entry is then 
ROKAP at the tangent point. If this again 
followed by zeros, linear interpolation is used 
to the next non-zero entry to yield ROKAP along 
a conical afterbody.) 


GROUP 5 STAGNATION DATA 
CARD 1 FORMAT (E10. 4) 


FIELD 1 (Columns 1-10) PTET. Local stagnation pressure 
in atmospheres. 


GROUP 6 FIRST GUESS INFORMATION 
CARD 1 FORMAT (5 (2XE14.7)) 


PRED 1 ,2,.-2.5 per card(16 columns per field) F(t.) 
First guesses for stream function (F(1,J)), 
velocity ratio (F(2,J)), shear function 
(F(3,J)), and derivative of shear function 
(F(4,J)) across the boundary layer. 


GROUP 7 STREAMWISE DISTRIBUTIONS FOR EDGE CONDITIONS 
CARD 1 FORMAT (5 (2XE14.7) ) 


PrELD 1,25...<5 per Card(l6 columns per field) PRE(IS). 
Ratio of local static to stagnation pressure 
for the mean flow condition. When KR(16) is 
non-zero PRE is input in the form of negative 
pressure coefficients referred to PTET. 


CARD 2 FORMAT(Blank) This card is inserted to maintain 
data continuity. In more complex formulations 
this Card corresponds to a control flag for 
updating the pressure profile during the com- 
PueEaeLon. 


Si 


GROUP 8 STREAMWISE DISTRIBUTIONS FOR WALL INPUT CONDITIONS 
CARD 1 FORMAT(E10.4) 


FIELD 1 (columns 1-10) T. Wall temperature in degrees 
R (for the present this is the free stream 
static temperature). 


CARD 2 FORMAT(8E10.4) USED ONLY FOR KR(9)=0 


FLELD 1L,2,...8,per card(l10 columns per field) FW(Uisoe 
Wall stream function (negative for mass 
addi tien). 


CARD 3 FORMAT(Blank) USED ONLY IF CARD 2 IS USED. This 
card is inserted to maintain data continuity. 
In more complex formulations this card cor- 
responds to a control flag for updating the 
wall stream function profile during the compu- 
tation. 


CARD 4 FORMAT(8E10.4) USED ONLY FOR KR(9)=1 


FIELD 1,2,...8 per card(10 columns per field) RHOVW(@tSe 
Total mass flux at the wall (LB/SEC FT**2 or 
dimensionless for KR(8)=0 or 1 respectively). 
These values are positive for mass injection. 


CARD 5 FORMAT(Blank) USED ONLY IF CARD 4 IS USED. This 
card is inserted to maintain data continuity. 
In more complex formulations this card cor- 
responds to a control flag for updating the 
total mass flux profile during the computation. 


GROUP 9 FREE STREAM DATA 
CARD 1 FORMAT (3E10. 4) 


FIELD 1 (Columns 1-10) UMFS. Mean free stream velocity 
in FT/SEC. 


FIELD 2 (Columns 11-20) FREQ. Frequency of mean free 
stream oscillations in HZ. 


FIELD 3 (Columns 21-30) EPS. Ratio of the amplitude of 
the free stream oscillation to the local mean 
free stream velocity. 


LAST CARD FORMAT(Al) The purpose of this entry is to per- 
mit a test on whether or not a new case is to 
follow. In the event a case does not converge 
in the alloted number of iterations, any re- 
maining cards for that case are read and then 
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ignored until either a period or a comma is 
encountered in column 1. A comma in column 1 
indicates that another case is to follow, while 
a period in column 1 indicates that this is the 
last card in the input deck. 
The preceding input instructions should permit the user 
to fully exercise the UBLIMP program with a minimum of 
difficulty. For more detailed information on various rou- 


tines the reader is referred to reference 18, which is the 


user manual prepared for the BLIMP program. 
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